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ABSTRACT 


It  is  relatively  easy  to  generate,  by  digital  computer, 
large  numbers  of  seemingly  independent  random  numbers  with  a 
uniform  distribution  over  a  fixed  range,  say  ~  ^  ^ 

Methods  of  generating  gaussian,  or  normal,  random  numbers  generally 
are  based  on  either  non-linear  transformations  on  random  numbers 
from  a  uniform  population,  or  the  summing  of  enough  independent 
numbers  from  a  uniform  population  for  the  central  limit  theorem 
to  be  applicable.  In  the  first  case  a  time-consuming  evaluation 
of  a  complicated  function  is  involved.  The  second  method  is  also 
slow  because  a  large  number  of  uniform  random  variables  must  be 
generated  and  summed  for  each  normal  random  variable  obtained. 

This  note  discloses  a  method  based  on  the  central  limit  theorem, 
except  that  the  summing  of  N  uniform  random  variables  gives  N 
normal  random  variables.  The  approach  is  to  form  an  N  dimensional 
vector  whose  components  are  uniform  random  variables,  multiply  the 
vector  by  a  Hadamard  matrix,  and  use  the  resulting  components  as 
normal  random  variables.  It  can  be  shown  that  the  resulting  N 
components  have  a  uniform  density  inside  a  N-dimensional  hypercube 
aligned  with  diagonals  along  the  coordinate  axes.  However,  the 
one  dimensional  marginal  densities,  the  two  dimensional  marginal 
densities,  indeed  all  the  marginal  densities  tend  toward  the  normal 
density  as  N  gets  large.  Furthermore  the  components  are  uncorre¬ 
lated  and  have  equal  variance,  independent  of  N.  However,  some 
of  the  fourth  moments,  which  should  be  zero  for  independent 
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normal  random  variables,  are  not  zero  for  our  derived  set 
(although  these  moments  do  approach  zero  as  N  becomes  large) . 
These  moments  can  be  made  zero,  however,  by  randomly  changing, 
or  not  changing,  the  sign  of  each  component. 

The  method  proposed  is  very  fast  because  the  principal 
step,  Hadamard  matrix  multiplication,  requires  only  N  log2  N 
additions  to  produce  N  components. 


Accepted  for  the  Air  Force 

Franklin  C.  Hudson 

Chief,  Lincoln  Laboratory  Office 
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A  NEW  METHOD  OF  GENERATING  GAUSSIAN  RANDOM 
VARIABLES  BY  COMPUTER 

Many  scientific  computations  require  the  generation  by  a 
computer  of  a  large  number  of  seemingly  random  numbers^  with  a 
multidimensional  normal  (gaussian)  probability  density  function. 
However,  although  there  are  efficient  algorithms  available  for 
generating  uniformly  distributed  random  numbers,  the  algorithms 
which  generate  normal  random  numbers  are  generally  much  slower. 

One  typical  algorithm  begins  by  generating  a  moderate  number  of 
uniform  random  variables  (say  N)  and  forming  the  sum,  which,  for 
large  N,  is  approximately  normal,  according  to  the  central  limit 
theorem.  This  method  is  slow  because  many  uniformly  distributed 
random  variables  must  be  developed  for  each  resulting  normal 
random  variable,  and  because  many  additions  must  be  performed 
for  each  normal  random  variable,  as  well.  Another  typical  method, 
more  attractive  because  it  is  less  dependent  on  approximation, 
begins  with  a  single  uniform  random  variable  and  performs  a  non¬ 
linear  transformation  of  the  random  variable  to  form  a  new  random 
variable  with  a  normal  density.  Unfortunately,  the  non-linear 
function  which  must  be  evaluated  is  complicated  and  cannot  be 
computed  quickly  on  most  digital  computers. 

This  note  describes  a  variation  of  the  central  limit  theorem 
approach,  in  which  N  uniform  random  variables  (r.v.)  are  trans¬ 
formed  into  N  approximately  normal  r.v.  by  a  Hadamard  transformation. 

T"!  The  term  "pseudorandom”  is  often  used  to  describe  the  determin¬ 
istically  generated,  but  seemingly  random  numbers. 
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This  operation  takes  y  N  log2  N  additions  and  y  N  log2  N  subtractions, 
which  is  to  say,  log2  N  additions  or  subtractions  per  normal  r.v, 
obtained,  and,  of  course,  it  is  only  necessary  to  generate  one 
uniform  r.v.  for  each  normal  r.v.  obtained.  The  plan  of  the 
paper  is  to  first  describe  the  ideal  uniform  and  normal  multivariate 
densities,  second  to  describe  the  Hadamard  transformation  used, 
third  to  derive  the  properties  of  the  transformed  variables,  and 
fourth,  to  consider  the  properties  of  the  transformed  variables 
after  they  have  been  subjected  to  random  sign  changes. 


I .  Properties  of  the  Ideal  Uniform  and  Normal  Densities 


We  define  the  ideal  uniform  r.v.,  X^,  as  one  whose  probability 


density  is 


Px(^n)  = 


1 

0 


Xj  < 


n 


(1) 


For  several  such  r.v.s  Xq,Xj^,  .  .  .  to  be  independent  implies 

that  the  joint  density  is 

1  if  all  |X^1  <  ^ 


(Xj-),Xi  ,X^  >  •  •  • 


XX. . .x^^O’^1’^2 


0  if  any  |x^l  >  j 


(2) 


On  an  N  dimensional  space,  the  joint  density  is  unity  on  the 
inside  and  zero  on  the  outside  of  a  hypercube  centered  on  the 
origin  with  coordinate  axes  passing  from  the  center  of  each  face 
to  the  center  of  the  opposite  face. 
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The  expected  value,  E(X^)  “ 

all  the  other  moments  can  be  found. 

moments  of  X  , 
n  ’ 


is  zero  and,  by  integration. 
For  the  case  of  even  order 


5?P  =  ^ 

^  2^P(2p+l) 


and  for  odd  order  moments 


(3) 


=  0 


(4) 


Since  the  mean  is  zero,  the  concept  of  central  moments  is  super¬ 
fluous.  For  joint  moments,  we  appeal  to  independence,  which  says 
that  the  moment  of  a  product  is  the  product  of  the  moments. 
Therefore 


=  (X^XXbXX^)  (5) 

and,  of  course,  these  joint  moments  are  zero  except  when  m,  n,  I 
are  all  even. 

It  is  worth  taking  note  of  the  second  moment,  =  1/12. 

The  normal  probability  density  of  an  r.v.  Y^,  with  zero  mean. 


is  given  by 
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where  a  is  the  variance.  We  shall  be  comparing  several  densities 
in  this  note  and  we  shall  force  the  variances  to  all  be  equal  to 
1/12.  Therefore  =  1/12. 

The  multivariate  density  of  several  such  independent  variables, 
•  •  •  >^N-l’ 


/2ti 


1  -ao+Yi+...+YN-l)/2o^ 


TT  a 


(7) 


The  odd  order  one  dimensional  moments  of  Y  are 

n 


=  0 
n 


(8) 


while  the  even  order  moments  may  be  shown  to  be 


p  ^  l-3»5* • • (2p-l)  _  l-3»5* • • (2p-l) 


(l/c^) 


I2P 


(9) 


Again,  by  independence,  the  joint  moments  are  equal  to  the  product 
of  the  individual  moments,  so  that  the  moment 


a'  ^  b^  ^  c • 


(10) 


is  zero  unless  m,  n,  I  are  all  even. 
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The  normal  density  is  never  zero,  and  is,  in  fact,  radially 
symmetrical,  i.e.,  Pyyy  y ^^O’^l’ *  *  * ’^N-1^  only  a  function  of 
the  distance  of  Yq,Y2^,  .  .  .  from  the  origin  in  N-space.  This 

is  also  a  true  statement  for  any  marginal  density.  The  same 
statement  is  not  true  for  the  multivariate  uniform  density,  which 
takes  on  different  probabilities  in  the  directions  of  the  vertices 
of  the  hypercube  than  along  the  direction  of  the  faces  of  the  cube. 

II.  The  Hadamard  Matrix 

By  definition,  a  Hadamard  matrix  is  an  orthogonal  matrix 
whose  elements  are  all  either  +1  or  -1.  Hadamard  matrices  do  not 
exist  for  all  dimensionalities;  for  example  there  is  no  3  x  3 
Hadamard  matrix.  It  is  relatively  easy  to  construct  an  N  X  N 
Hadamard  matrix,  however,  if  N  is  a  power  of  two.  We  shall  content 
ourselves  with  a  particular  Hadamard  matrix  for  each  power  of 
two,  for  our  purpose  is  not  to  investigate  Hadamard  matrices  but 
to  apply  them.  Letting  N  =  2  ,  and  numbering  our  dimensions  in 
the  binary  system 

i  =  i^  +  2ii  +  4io  +  . . .  +  2^"^i  .  0  i  s  N-1 

U  i  z  c-i 

(11) 

j  =  jo  +  2Jl  +  ^^2  +  •••  ■*■  ^‘'■’■jc-l  0  s  j  i  N-1 

the  elements  of  our  particular  Hadamard  matrix,  H,  are  given  by 
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hj^  .  -  (-1)’'°^°  (-1)’'’-^^  (-1)''^^^ 


(-1) 


^c-l^c-1 


(12) 


We  shall  show  in  an  appendix  that  multiplication  of  an  N  component 
vector  by  H  can  be  accomplished  with  N  log2  N  additions  and 
subtractions . 

It  is  clear  that  h.  .  =  h.  .  so  that  H  is  a  symmetric  matrix. 

^ » J  J  » ^ 

ttl  ttl 

Also,  the  zero^  row  and  the  zero^  column  are  composed  of  +l's 
only.  The  product  of  the  elements  in  the  i^  row  with  the  row 

give  the  elements  of  another  row  of  the  Hadamard  matrix. 


k,j 


=  h 


where  i  is  the  row  whose  number  is 


(13) 


X  =  (iQ  ©  k^)  +  2(i^  ©  k^^)  +  ... 

+  e  k^.i)  (14) 

where  ©  means  "exclusive  or".  We  shall  say  as  a  shorthand, 

jJ  =  i  ©  k 

in  place  of  (14).  By  symmetry,  the  product  of  columns  is  also  a 
column. 
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A  consequence  of  (13)  is  the  orthogonality  of  the  Hadamard 


matrix.  Strictly  speaking  the  matrix 


1 

7n 


H 


is  not  only  its  own  transpose  but  its  own  inverse,  and  it  is  the 
matrix  which  we  shall  use  in  the  transformation  of  uniform  r.v. 
to  normal  r.v. 

We  give  as  an  example  the  4x4  Hadamard  matrix 


H 


1 

1 

1 

1 


1  1 

-1  1 

1  -1 

-1  -1 


1 

-1 


-1 

1 


III.  Hadamard  Generation  of  Almost  Normal  r.v. 


We  begin  with  a  set  of  uniformly  distributed  r.v.  assumed 
independent.  We  form  a  new  set  of  r.v.  from  the  uniform  set  by 
the  operation 


W  —  -7^ 
m  /N 


N-1 
S  h, 
n=0 


X 


n,m  n 


m  =  0,1, 


.N-1 


(15) 


and  we  shall  be  interested  in  the  properties  of  the  set  of  r.v. 

W^.  For  any  single  random  variable,  it  is  easy  to  see  that  P^(Wj^) 
is  the  density  of  a  sum  of  independent  r.v.  with  zero  mean  and 
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equal  variance.  It  is  known  from  the  central  limit  theorem  that 
such  a  sum  approaches  a  gaussian  r.v.  in  the  limit  of  large  N. 


The  approximation  is  good  in  the  region  of  near  zero,  bad  in 
the  tails.  It  is,  of  course,  clear  that 


(16) 


so  that  Is  zero  outside  this  range.  However,  for  large  N, 

^  /N  is  many  standard  deviations  from  the  mean,  so  that  the 
difference  can  be  tolerated.  Indeed,  the  approximation  is  the 
same  as  for  the  common  method  of  generating  normal  r.v.  except 
that  we  hope  to  use  a  larger  value  of  N  and  so  obtain  a  good 
approximation . 


We  have  generated,  however,  N  different  r.v.  from  the  same 


N  uniform  r.v.  and  we  should  not  expect  the  to  be  independent. 

This  leads  us  to  investigate  the  multivariate  density 

^ww. .  .w^^O’^1’ •  •  • ’^N-1^  • 

(wqj^I,  .  .  .  =  ^xx.  .  .x^^O’^1’ •  " ’^N-1^/ 


ww. 


.w 


where 


D 


3W^ 

^^o 


£ 

N-1 


BW. 


N-1, . . . , 

0 


BW 

TK 


N-l 

N-1 


P 

ww. . .w 


(Wq,W3_, 


>Wn-i) 


p 

XX ...  X 


•  • 


' ’^N-l 


) 


(17) 
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since  the  magnitude  of  the  determinant  of  ^  H  is  unity.  Of 
course,  in  (17),  P  (X^,X,  ,  .  .  .  ,X,,  ,)  is  to  be  evaluated  at 


"  ^0  ■ 
^1 

1  ,T 

■  ^0  ■ 
Wi 

• 

1 

X 

> 

t  _ 

f 

•  • 

7* 

mA 

1 _ 

but  this  is  only  a  change  of  coordinate  system.  In  fact,  there¬ 
fore,  the  N-dimensional  joint  density  P  is  not  normal  but 

uniform.  The  density  is  one  inside  and  zero  outside  a  hypercube, 
This  hypercube,  however,  is  rotated.  Whereas  for  the  variables 
X^,  the  hypercube  was  situated  with  the  coordinate  axes  passing 
through  its  faces,  for  the  W^'s  we  shall  see  that  the  coordinate 
axes  pass  through  vertices  of  the  hypercube.  In  fact,  the 
coordinates  of  some  of  the  vertices  of  the  hypercube  surrounding 

^xx. .  •  •  • ’\-l^ 


Xj  =  j 


h .  . 


for  any  i,  and  such  a  vertex  becomes  rotated  into  the  position 
with  coordinates 


W. 

J 


0 


J  =  1 

j  ^  i 
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tin 

which  lies  on  the  coordinate  axis  of  the  i^  coordinate  in  W-space. 
The  vertex  at  -X^  rotates  to  the  opposite  position,  thus  proving 
that  the  coordinate  axes  of  W-space  pass  through  diagonals  of  the 
hypercub^^ . 

A  hypercube  in  N-dimensions  has  2^  vertices.  There  are, 

however,  only  N  coordinate  axes,  so  only  2N  vertices  can  lie  on 

the  coordinate  axes.  The  remaining  2^  -  2N  vertices  of  the 

hypercube  do  not  map  onto  the  coordinate  axes.  At  least  some  of 

the  vertices  are  found  in  the  "interior"  of  hyperquadrants,  but 

not  every  hyperquadrant  of  the  N-space  can  be  lucky  enough  to  have 

a  vertex  within  it,  for  there  are  2^  hyperquadrants  and  at  most 

2^  -  2N  vertices  to  go  around.  Therefore,  we  learn  that 

w^^0*^l*  *  *  *  *^N-1^  only  radially  as5rmmetrical  but  is 

not  S5rmmetrical  in  each  hyperquadrant  of  W-space.  This  shows  up 

when  we  compute  the  moments  of  P  .  However,  the  most  important 

ww . . .w  ’ 

moment  of  P  ,  namely  W  W  ,  is  zero,  showing  that  the  variables 

ww .  .  .w*  m  n*  *  ^ 

W  are  at  least  uncorrelated.  We  shall  now  compute  all  the  first 
m 

and  second  moments  of  P,„,  (Wn»Wi  ,  .  .  .  ,W..  ,). 

WW...W'  U’  i  ’  N-i 

The  first  moments  are  all  equal  and  are  all  zero. 


W  =  -TT? 
m  /N 


E  h 

r 

n 


E  h 

I 

n 


\  =  0 


(18) 
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because  X  =0.  The  second  moments  are 
n 


N-1  N-1 


WW,  =  Eh  h 

a  b  N  _  rv  r\ 
m=0  n=0 


,a  n,b  Vn 


(19) 


but  is  zero  unless  m  =  n,  in  which  case  it  is  a  .  Thus  (19) 

becomes 


_ _  1  N-1  ^ 

W  W,  =  ^  E  h  /  ffluN  o 

a  b  N  m,  (a®b) 

m=0 


(20) 


th 


By  definition  of  a  ©  b,  it  is  zero  if  a  =  b,  and  the  zero  column 
of  a  Hadamard  matrix  consists  of  N  ones.  Therefore 


^  2 
W  =  a 
a 


(21) 


as  expected.  But  if  a  7^  b,  a  ©  b  is  not  zero  and  sirce  all  the 
other  columns  of  H  contain  as  many  -I's  as  I's,  the  sum 


E  h 
m 


m,  (a®b) 


IS  zero,  giving 


W  W, 
a  b 


=  0 


as  stated  above. 


(22) 
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It  is  not  difficult  to  argue  that  all  the  odd  order  moments 

of  „  are  zero.  We  focus  now  on  fourth  moments.  By  definition 

ww  •  .  .  w 


W  W,  W  W  ,  =  -TT 
a  b  c  d 


h.  h.  ,h 
j  ,b 


N-1 

N-1 

N-1 

N-1 

Z 

Z 

Z 

Z 

i=0 

j=o 

O 

II 

£=0 

1  ,  X.X.X,  X, 

k,c  jJ,d  1  j  k 


(23) 


To  evaluate  the  fourth  moments  we  note  that  X.X.X,  X.  =  0  for 

1  j  k 

most  combinations  of  i,j,k,;i  with  the  exceptions  listed  below: 

1 


1. 

2. 

3. 

4. 


O’ 


a  = 


WT 


i=j  ,  V.=  SLy  i/^k  for  which  X^X^  = 

— T  L  1 

i=k,  j  =  Ji,  t/i  for  which  =  o  = 

i=;i,  j=k,  iT^j  for  which  X?X?  = 


i=j=k=£ 


“ZT 

for  which  X. 


1 


These  are  all  disjoint  cases.  Therefore 


W 


1  ,  1 


2/n-I 


/N-1 


aVc^d  ==  ^i,a0tij^fQ  \,cedl 

Vk^i 

^  .  2/n-1  /n-1  \ 

^  K  *’i.a0c(  jfo  \.bOd 


1  / 1 


2/n-1 


'n-1 


U=0  ''i.awL.o 


N-1 


^  ^i,a0b9=0d 


(24) 
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Equation  (24)  looks  very  complicated  but  it  may  be  straightforwardly 

-T 

evaluated  for  any  special  case.  Thus,  for  W  we  have  a=b=c=d  and 

3. 

we  f ind 


^  (N^-N)  +  (^)  (N^-N)  + 

a  iz  iz 

(^)  (N^-N)  + 

N  N 

^  ^  TM  ~  ■  TM  <^25) 

Similarly,  if  we  let  a=c ,  b=d,  aj^b  we  get 

“a?  '  ^  (N^-N)  +  ^  (^)N 

^  '  IW  +  ISON  =  (26) 

These  results  approach  the  gauss ian  independent  results 

for  large  N.  But  consider  the  curious  case  of  a,b,c,d  all  unequal 

but  with  a0b0c0d  =  0.  For  this 


^a^b^c^d  =  m 


(27) 
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An  example  of  such  a  fourth  moment  is  This 

moment  also  approaches  the  expected  result  for  independent  normal 
r.v.  as  N  However,  it  is  somewhat  disturbing  that  it  is 

non-zero.  We  can,  after  all,  convince  ourselves  that  if  a 
certain  moment  is  within,  say,  27o  of  the  expected  result  for 
normal  independent  r.v.  that  this  is  somehow  good  enough.  But 
what  can  we  say  if  the  expected  moment  should  be  zero?  It  is 
not  clear  how  important  this  is,  and  indeed  it  probably  depends 
on  the  application  for  which  the  r.v.  are  desired. 

Before  considering  the  effect  of  a  random  sign  change,  it 
is  well  to  summarize  what  we  have  shown  so  far.  We  have  proposed 
a  method  by  which  uniformly  distributed  zero-mean  random  noise 
can  be  converted  to  approximately  normal  zero  mean  random  noise. 
The  one  dimensional  marginal  densities  are  sums  of  N  independent 
uniform  r.v.  and  thus  tend  nicely  toward  normal  for  large  N. 

The  r.v.  are  also  uncorrelated  (in  the  sense  that  the  expected 
value  of  a  product  is  zero).  Furthermore,  based  on  the  moments 
we  have  computed,  the  difference  between  a  given  joint  moment 
and  the  same  joint  moment  for  independent  gaussian  noise  seems 
to  go  to  zero  as  1/N.  N  can  reasonably  be  made  quite  large  since 
the  amount  of  computation  is  proportional  to  log2  N  per  normal 
r.v.  generated. 
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IV.  HADAMARD  GENERATION  OF  MORE  NEARLY  INDEPENDENT  RANDOM 

VARIABLES 


Let  us  now  suppose  that  we  have  generated  N  random  variables 

W  .  We  have  seen  that  some  of  the  moments  W  W,  W  W,  ,  which 
m  abed’ 

would  be  zero  for  an  independent  gauss ian  set,  are  non-zero.  We 
can  form  a  set  of  N  new  random  variables  from  the  W's  by  changing 
the  sign  of  each  W,  or  not  changing  it,  at  random.  That  is, 
consider  an  easily  obtained  set  of  r.v.,  r^  ,  independent  of  each 
other  and  of  W^,  with  equal  likelihood  of  being  +1  or  -1.  For 
these  r.v.  the  moments  are  all  either  zero  or  one.  Even  order 
moments  are  equal  to  1, 


odd  order  moments  are  equal  to  zero. 


2p+  1 

r  ^ 


=  0 


5 


and,  by  independence,  the  moment  of  the  product  of  different 
is  the  product  of  the  moments.  Now  consider  the  set  of  random 
variables 


V  »  r  W 
m  mm 


(28) 
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We  shall  compute  some  of  the  moments  of  the  joint  density 


^vv  ^l’*'''^N“l^  show  that  this  density  is  nearly 

gauss ian,  for  large  N.  First  consider  the  moments  which 


violated  the  gaussian  assumption  for  the  W's.  The 
equal  to 


(WaVc«d) 


and  unless  a,  b,  c,  d  are  equal  at  least  in  pairs,  ^^^^^c^d 
be  zero.  A  similar  argument  will  show  that  all  joint  moments  of 

P  which  would  be  zero  for  independent,  zero  mean  gaussian 

noise  are  indeed  zero.  However,  all  other  moments  will  be 

— I  2 

unaffected  by  the  conversion  of  W  to  V  .  Thus  V  =  a  , 

m  m  a  ’ 

V  V,  =  0,  and 
a  b  ’ 


„  2„  2  „  2„  2  4  ,  a 

''a  ''b  *  "a  "b  ■  "  +  T5N 


(29) 


We  next  consider  the 


order  moments 


Discarding  the  trivially  zero  moments,  we  have  for  the  other 


a  D  c  d  e  f 


1 


E  E  L  E  L  L 
i  j  k  m  n 


h.  h.  , 
i.a  j,b 


h,  j 
C  1,0. 


h  h  .  X.X.X,  X.X  X 
m,e  n,f  ijkJimn 


We  must  again  consider  moments  X.X.XiX.X  X  ,  most  of  which  are 

ijKj6mn’ 

zero.  The  exceptions  are  listed  below  with  the  equal  subscripts 
grouped  and  the  ungrouped  subscripts  unequal.  Thus 


Disjoint  Cases 


Moment 


i  j  k  A  m  n 


1 

448 


i 

j 

k 

1 

m 

n 

i 

j 

k 

m 

1 

n 

i 

j 

k 

n 

1 

m 

i 

j 

1 

m 

k 

n 

i 

j 

1 

n 

k 

m 

1 

i 

k 

1 

n 

j 

n 

i 

k 

1 

n 

j 

m 

j 

k 

1 

m 

i 

n 

j 

k 

i 

n 

i 

m 

k 

1 

m 

n 

i 

j 

j 

1 

m 

n 

X 

k 

j 

k 

m 

n 

i 

a 

i 

k 

m 

n 

j 

SL 

i 

1 

m 

n 

J 

k 

i 

J 

m 

n 

i 

i 
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Disjoint  Cases 


Moment 


i  j 
i  j 
i  j 

i  k 
i  k 
i  k 
i 

i  I 
i  I 
i  m 
i  m 
i  m 
i  n 
i  n 
i  n 


k  I, 
k  m 
k  n 

j  ^ 
j  m 

j  n 
j  k 
j  m 
j  n 
j  k 

j  ^ 
j  n 
j  k 

j  ^ 
j  m 


m  n 
n  -t 
m  I, 
m  n 
I,  n 
I,  m 
m  n 
n  k 
m  k 
I,  n 
k  n 
k  -t 
I,  m 
k  m 
k  -t 


1 

1728 


t  h 

This  table  allows  us  to  compute  all  the  6*"  order  moments 
for  special  cases.  For  example,  if  we  were  to  compute  V  ^  ,  we 

oL 

would  have  as  part  of  the  computation  a  sum 


1 


N-1 

L 

i=0 


hi, 


(a@a@a@a9a 


7^ 

0  a)  ^i 


1  1 
548 


and  fifteen  sums  similar  to 
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^  N-1 

0  a  ©  a  0  a) 


N-1 
""n  h 

m=0  m, 
mif^i 


a 


0  a 


1 

N-^ 

■  9^ 

fifteen  sums 

similar  to 

,  N-1 

N-1 

N-1 

~T  ^ 
N-^  i=0 

^i,  (a  0  a) 

L  h, 
k=0 

(a  ®  m=0 

k?^i 

mf^i 

mf^k 

2  2  2 


1 

1728 


t  h 

Gathering  all  the  terms  together,  we  get,  for  the  6*"  moment 


15 

1728 


N(N-l)(N-2)  ,  15  /N(N-l),.  1  ,  N  >. 

- ^3 - +  960  ^  448 


V 

a 


15ct°  +  terms 


in  ^  and 


(30) 


— 5 — 2^ 

Evaluating  V  V,  we  set,  for  example,  a  =  c  =  d  =  e,  b  =  f, 

3.  D 

a  /  b.  Therefore,  the  sxim 
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is  the  same  as  before,  but  many  of  the  other  thirty  sums  are 
zero,  leaving  only  the  cases 


i 

j 

1 

n 

k 

m 

i 

k 

j 

n 

1 

m 

i 

k 

1 

m 

j 

n 

i 

1 

j 

n 

m 

k 

j 

k 

1 

n 

i 

m 

and 

i 

m 

j 

n 

k 

1 

j 

1 

m 

n 

i 

k 

j 

k 

m 

n 

i 

1 

i 

j 

m 

n 

k 

1 

We  quickly  count 


3 

TT^ 


(  N(N-l)(N-2)^ 

IT 


960 


,N(N-1).  ,  1  ,N  >, 


+  terms 


•  1  1 
in  ^  and  — »y 
N  ^ 


(31) 


— 5 — 5 — 5" 

Similarly  for  V  V,  we  obtain 

a  D  c 


— 2 — 5 — 5" 

V  V,  V 
a  b  c 


— 2 — 2 — 2" 

a  b  c 


_  .  N(N-l)(N-2).  3  (N(N-l).  1  ,N  ) 


+  terms  in  h  and  i 


(32) 


These  results  also  approach  the  expected  results  for 

6  6  6 

normal  r.v.  with  an  error  like  1/N,  i.e.,  15a  ,  3a°,  and  a°. 

At  this  point,  we  comment  that  the  behavior,  for  large  N, 
of  any  moment  is  a  matter  of  combinatorial  analysis.  For  example, 
the  moment 


^  =  1 *3 *5 • • •  (2p-l) +  terms  in  ^  ,  etc.  (33) 

can  be  derived  by  a  counting  process.  For  more  complicated 
moments,  the  analysis  is  so  complicated  that  we  have  not  under¬ 
taken  it. 

We  previously  described  the  N  dimensional  distribution  of 

the  W's  in  terms  of  a  hypercube  which  we  argued  was  asjnranetically 

N 

distributed  among  the  2  hyperquadrants  in  N  space.  We  cannot 
make  a  simple  statement  about  ^  ^^0’  ^l’‘''^N-l^  because 

it  is  not  uniformly  distributed,  but  rather  has  some  density 
which  is  a  continuous  function  of  spacial  coordinates.  However, 
it  is  clearly  symmetrically  distributed  in  the  hyperquadrants, 
by  construction. 

Programming  Considerations 

The  method  proposed  requires  at  least  the  following  two  or 
three  steps: 
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1.  Compute  N  uniform  r.v.  Xq,  Suppose  that  each 

r.v.  requires  a  time  to  generate.  Then  part  1  will  require 
a  time  NT^. 

2.  Compute  the  "Hadamard  transform”  of  the  array  of  N  numbers 

from  part  1.  This  can  be  done  in-place  (i.e.,  without  any  extra 

memory)  with  N  log2N  additions  or  subtractions.  Assuming  that 

an  addition  and  a  subtraction  take  the  same  amount  of  time  T  , 

a’ 

part  2  will  require  a  time  N(log^N)T  .  If  it  is  desired  to 

Z  o. 

adjust  the  variance  of  the  resulting  W^,  we  should  add  to  this 
time  an  additional  time  N  T  ,  where  T  is  the  time  for  a  multiply 
(or  a  scale  if  that  is  acceptable) . 


3.  If  the  improved  approximation  to  independent  gaussian  noise 

is  desired,  generate  N  additional  random  variables  r^  which  can 

be  used  to  randomly  change  or  not  change  the  sign  of  the  r.v. 

Wj  generated  in  part  2.  This  should  take  at  most  a  time 

N(T  +  T  )  assuming  the  r.  can  be  generated  as  quickly  as  a 
u  a  j 

uniform  r.v.  and  sign  changing  can  be  accomplished  as  quickly 
as  addition  or  subtraction. 

The  total  time,  therefore,  is  like 


’'total  '  +  log2(2N)  T^) 

but  this  is  the  time  for  generating  N  r.v.  whereas  the  relevant 
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time  would  be  the  time  per  random  variable.  This  latter  time 
is 


Tg  =  +  <1°82”  +  ' 

If  many  r.v.  are  needed,  as  will  ixsually  be  the  case,  they 
should  be  generated  N  at  a  time  by  a  subroutine.  This  sub¬ 
routine  could,  of  course,  supply  them  one  at  a  time  to  a  user 
program,  recomputing  a  block  of  N  when  its  supply  was  exhausted. 
Thus,  unless  fewer  than  N  r.v.  were  required,  the  limitation  on 
computing  the  r.v.  N  at  a  time  should  not  be  important. 

Conclusions 

Random  variables  with  an  almost  normal  distribution  can  be 
computed  by  first  obtaining  uniform  random  variables,  using 
standard  computer  algorithms,  and  then  performing  a  linear  trans- 
formation  on  a  set  of  N  =  2  of  the  uniform  random  variables  to 
obtain  N  random  variables  whose  joint  distribution  is  not 
gaussian  but  uniform  in  a  hypercubic  region.  Nevertheless,  many 
of  the  marginal  densities  are  nearly  gaussian,  in  accord  with 
the  central  limit  theorem,  and  the  r.v.  are  uncorrelated.  In 
the  limit  of  large  N,  all  the  moments  which  have  been  determined 
seem  to  approach  the  moments  expected  for  gaussian  independent 
r.v.  However,  some  of  the  joint  moments  which  would  be  zero  for 
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independent  gaussian  r.v.  are  finite  for  the  derived  set.  By 
randomly  changing  or  not  changing  the  sign  of  the  linearly  trans¬ 
formed  r.v. ,  a  new  set  of  r.v.  is  obtained  for  which  all  the 
moments,  including  joint  moments,  show  considerable  similarity 
to  the  moments  expected  from  the  gaussian  independent  case,  for 
large  N.  The  procedure  has  the  advantage  of  being  faster  than 
other  decent  methods  of  generating  gaussian  r.v.  on  a  digital 
computer,  with  relatively  good  accuracy. 
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APPENDIX 


NUMBER  OF  OPERATIONS  IN  HADAMARD  MATRIX  MULTIPLICATION 

If  we  consider  the  multiplication  of  an  arbitrary  vector  by 

the  matrix  H,  with  elements  h.  .  given  by  (12),  an  abvious 

^  I  J 

proceedure  is  to  partition  the  matrix  into  four  quadrants, 


H  = 


c .  . 


b.  . 


=  h.  .  -N 

“ij 

=  h.  .  .  N 

+  J 

0  ^  i  ^  I  - 

=  h.  ,  N  . 

1  +  j,J 

0  S  J  S  1  - 

'‘i.J 

=  h.  ,  N  .  ,  N 

1  +  j.J  +  5-  ^ 

1  that 

for 

a .  .  • 

^0^  0  ^1^  1 
=  (-1)  ^  ^  (-1)  ^  ^  ... 

(.l)ic-2jc-2 

(-1) 


-l^c-l 


i^_^  jc-1  ~  ^  entire  quadrant.  Therefore  a^  ^  is 


actually  an  ^  x  y  Hadamard  matrix.  Similarly 


a  •  •  “*  b  •  •  c  •  •  “*  "d  •  • 

3-,j  i,J 
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It  follows  that  if  the  vector  to  be  multiplied  by  H  in  a 
column  vector  [X^l  we  may  form  two  shorter  vectors  with  components 


I 

^i 


It 

^i 


X.  +  X.  ,  N 

1  ±  ’T 


X.  -  X.  ,  N 
1  ^  7 


0  ^  i  ^ 


0  i  i  i 


N 


N 


-  1 


-  1 


N 

and  we  shall  find  that  the  first  components  of  the  product  are 
given  by 

N 


Z 

i=0 


a .  . 


N 


X. 


and  the  second  components  are  given  by 
N 


E  a .  . 
i=0 


It 

X.. 


N  ..  N 


each  of  which  are  ■j"  ^  'J'  Hadamard  matrix  multiplications.  This 
reduction  cost  us  N/2  additions  and  N/2  subtractions.  Of  course 


a.  .  may  be  similarly  partitioned  and  this  may  be  done  repeatedly. 
^  >  J 

N  N 

Each  reduction  .costs  and  additional  ^  additions  and  ^  subtractions 
so  that  a  total  of  Me  =  N  log2  N  operations  is  required  in  all. 
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